查看原文
其他

基于 R 语言的科研论文绘图技巧详解(5)

庄闪闪 庄闪闪的R语言手册 2023-04-03

点击下方公众号,回复资料分享,收获惊喜

简介

在查阅文献的过程中,看到了几幅非常不错的出版图,今天就跟着小编一起学习下,他们是怎么使用 R 绘制出来的。

今天主要介绍 第五幅图(E) —— 带置信区间的拟合曲线图及带误差项的散点图。这个图在科研绘图中较为常用,例如:绘制产品失效时间的散点图。

前四幅图的详细代码介绍可见:基于 R 语言的科研论文绘图技巧详解(4)基于 R 语言的科研论文绘图技巧详解(3)基于 R 语言的科研论文绘图技巧详解(2)基于 R 语言的科研论文绘图技巧详解(1) 。最后一幅图会随后继续介绍,读者在学习过程中,可以将内部学到的知识点应用到自己的图形绘制中。推文已经将主要知识点进行罗列,更有利于读者学习和查阅。

那我们来看看,作者是怎么实现这个功能的吧,本文知识点较多,大家耐心学习,建议自己实践。对应代码、数据可在 GitHub - marco-meer/scifig_plot_examples_R: Scientific publication figure plotting examples with R[1] 中找到。

主要知识点

  • 绘制散点图、丝带形状图;
  • 绘制横向、纵向误差图;
  • 学会小技巧:展示轴外部的图形。

绘图

加载包

首先加载一些需要使用到的包。

library(ggplot2) # Grammar of graphics

设置主题

接下来,为了方便起见,作者在绘图前设置好了主题,并将该函数命名为 my_theme

这一部分在第一篇推文给出,代码将在文末中完整代码给出。

手动修改大部分面板,具体可以参考本篇文章[2]。或者观看我在 B 站发布的《R 语言可视化教程》,里面也有一些简单主题设置介绍。

导入数据

首先使用 read.csv() 导入三个基因相关数据,最后合并到 data_E 数据框中。前 6 行数据如下所示。

data_Ea = read.csv("./data_Ea.csv")
data_Eb = read.csv("./data_Eb.csv")
data_Ec = read.csv("./data_Ec.csv")
data_E = data.frame("gene"=c(rep("gene a",nrow(data_Ea)),
                             rep("gene b",nrow(data_Eb)),
                             rep("gene c",nrow(data_Ec))),
                    "t"=c(data_Ea$t,data_Eb$t,data_Ec$t),
                    "C"=c(data_Ea$C,data_Eb$C,data_Ec$C),
                    "pos_err"=c(data_Ea$err,data_Eb$pos_err,data_Ec$pos_err_C),
                    "neg_err"=c(data_Ea$err,data_Eb$neg_err,data_Ec$neg_err_C),
                    "pos_err_t"=c(rep(NA, nrow(data_Ea)),rep(NA, nrow(data_Eb)),data_Ec$pos_err_t),
                    "neg_err_t"=c(rep(NA, nrow(data_Ea)),rep(NA, nrow(data_Eb)),data_Ec$neg_err_t)
)

head(data_E)
#    gene  t    C pos_err neg_err pos_err_t neg_err_t
# 1 gene a 12 0.27    0.03    0.03        NA        NA
# 2 gene a 24 0.19    0.02    0.02        NA        NA
# 3 gene a 36 0.17    0.01    0.01        NA        NA
# 4 gene a 48 0.13    0.04    0.04        NA        NA
# 5 gene a 60 0.09    0.03    0.03        NA        NA
# 6 gene a 72 0.10    0.01    0.01        NA        NA

gene 是分类变量,t 表示时间将作为 x 轴。C 表示正常水平下的结果(小编对基因方面不是很懂,可能理解的不对)。pos_errneg_err 是正向与反向误差。

自定义函数

首先,定义三个函数,以便稍后拟合散点图。这里的参数具体怎么得到的,作者没有介绍,小编不是很了解。或许是使用了某种参数估计的方法吧。

f1 = function(t) 0.2*exp(-t/48)
f2 = function(t) 0.3*exp(-t/60)
f3 = function(t) 0.4*exp(-t/72)
t = seq(0,96,1)

# 构造 f2 函数的数据
ribbon = data.frame(
  "f2" = 0.3*exp(-t/60),
  "t" = t
)

绘制图形

定义图形形状分别为实心正方形、圆形、三角形 manual_pch =c(15,16,17) 。之后简单绘制下散点图geom_point()和 丝带形状的图 geom_ribbon()

注意:这里的 factor() 将 gene 变成离散形式。由于 geom_ribbon() 中使用的数据不同,所以要重新指定下新数据。

manual_pch =c(15,16,17) 定义图形形状
ggplot(data=data_E) +
  geom_point(aes(x=t,y=C,pch=factor(gene)),size=2) +
  geom_ribbon(data=ribbon, aes(x=t,ymin=0.85*f2,ymax=1.15*f2),
              fill="black",alpha=0.1

在这个基础上使用 stat_function() 添加三个函数,两条虚线,一条实线。

stat_function(fun=f1, geom="line",linetype="dashed") +
  stat_function(fun=f2, geom="line") +
  stat_function(fun=f3, geom="line",linetype="dotted")

使用 geom_errorbar()geom_errorbarh() 添加误差棒(纵向与横向)。注意内部参数不同:(ymin,ymax) 和 (xmin,xmax)。

geom_errorbar(aes(x=t,ymin=C-neg_err, ymax=C+pos_err), width=2) +
geom_errorbarh(aes(y=C,xmin=t-neg_err_t,xmax=t+pos_err_t))

最后就是主题的设置了,包括运用了自定义的主题 my_theme()

scale_shape_manual(values=manual_pch) +
  my_theme() + theme(legend.title=element_blank())+
  theme(legend.position = c(0.9,0.95)) +
  scale_x_continuous(expand = c(0, 0), 
                     breaks = c(seq(0,96,12)),
                     limits = c(0,96)
  ) +
  scale_y_continuous(expand = c(0, 0),
                     breaks = c(seq(0,0.4,0.05)),
                     limits = c(0,0.4)
  ) +
  theme(
    panel.grid.major = element_line("gray95", size = 0.1),
    # putting label closer to axis bc exponent makes it bigger 
    axis.title.y = element_text(margin = margin(r = -9)) 
  ) +
  xlab("time (h)") +
  ylab(expression(paste("concentration (mmol",~cm^-3,")"))) 
 

但是,最右边几个点并没有显示清楚,这里作者运用了一个小技巧,将这些点清晰的呈现出来。来看看这个代码的效果:

 coord_cartesian(clip = "off"# to allow for plotting outside axes

完整代码

# Panel E ----
library(ggplot2) 
base_size = 12 
my_theme <-  function() {
  theme(
    aspect.ratio = 1,
    axis.line =element_line(colour = "black"),  
    
    # shift axis text closer to axis bc ticks are facing inwards
    axis.text.x = element_text(size = base_size*0.8, color = "black"
                               lineheight = 0.9,
                               margin=unit(c(0.3,0.3,0.3,0.3), "cm")), 
    axis.text.y = element_text(size = base_size*0.8, color = "black"
                               lineheight = 0.9,
                               margin=unit(c(0.3,0.3,0.3,0.3), "cm")),  
    
    axis.ticks = element_line(color = "black", size  =  0.2),  
    axis.title.x = element_text(size = base_size, 
                                color = "black"
                                margin = margin(t = -5)), 
    # t (top), r (right), b (bottom), l (left)
    axis.title.y = element_text(size = base_size, 
                                color = "black", angle = 90,
                                margin = margin(r = -5)),  
    axis.ticks.length = unit(-0.3"lines"),  
    legend.background = element_rect(color = NA
                                     fill = NA), 
    legend.key = element_rect(color = "black",  
                              fill = "white"),  
    legend.key.size = unit(0.5"lines"),  
    legend.key.height =NULL,  
    legend.key.width = NULL,      
    legend.text = element_text(size = 0.6*base_size, 
                               color = "black"),  
    legend.title = element_text(size = 0.6*base_size, 
                                face = "bold"
                                hjust = 0
                                color = "black"),  
    legend.text.align = NULL,  
    legend.title.align = NULL,  
    legend.direction = "vertical",  
    legend.box = NULL
    panel.background = element_rect(fill = "white"
                                    color  =  NA),  
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    plot.title = element_text(size = base_size, 
                              color = "black"), 
  ) 
  
}

data_Ea = read.csv("./data_Ea.csv")
data_Eb = read.csv("./data_Eb.csv")
data_Ec = read.csv("./data_Ec.csv")
data_E = data.frame("gene"=c(rep("gene a",nrow(data_Ea)),
                             rep("gene b",nrow(data_Eb)),
                             rep("gene c",nrow(data_Ec))),
                    "t"=c(data_Ea$t,data_Eb$t,data_Ec$t),
                    "C"=c(data_Ea$C,data_Eb$C,data_Ec$C),
                    "pos_err"=c(data_Ea$err,data_Eb$pos_err,data_Ec$pos_err_C),
                    "neg_err"=c(data_Ea$err,data_Eb$neg_err,data_Ec$neg_err_C),
                    "pos_err_t"=c(rep(NA, nrow(data_Ea)),rep(NA, nrow(data_Eb)),data_Ec$pos_err_t),
                    "neg_err_t"=c(rep(NA, nrow(data_Ea)),rep(NA, nrow(data_Eb)),data_Ec$neg_err_t)
)

head(data_E)

f1 = function(t) 0.2*exp(-t/48)
f2 = function(t) 0.3*exp(-t/60)
f3 = function(t) 0.4*exp(-t/72)
t = seq(0,96,1)

ribbon = data.frame(
  "f2" = 0.3*exp(-t/60),
  "t" = t
)
head(ribbon)

manual_pch =c(15,16,17# available pch: type ?pch 

panel_E <- ggplot(data=data_E) +
  geom_point(aes(x=t,y=C,pch=factor(gene)),size=2) +
  geom_ribbon(data=ribbon, aes(x=t,ymin=0.85*f2,ymax=1.15*f2),
              fill="black",alpha=0.1) +
  stat_function(fun=f1, geom="line",linetype="dashed") +
  stat_function(fun=f2, geom="line") +
  stat_function(fun=f3, geom="line",linetype="dotted") +
  geom_errorbar(aes(x=t,ymin=C-neg_err, ymax=C+pos_err), 
                width=2) +
  geom_errorbarh(aes(y=C,xmin=t-neg_err_t,xmax=t+pos_err_t))+
  scale_shape_manual(values=manual_pch) +
  my_theme() + theme(legend.title=element_blank())+
  theme(legend.position = c(0.9,0.95)) +
  scale_x_continuous(expand = c(00), 
                     breaks = c(seq(0,96,12)),
                     limits = c(0,96)
  ) +
  scale_y_continuous(expand = c(00),
                     breaks = c(seq(0,0.4,0.05)),
                     limits = c(0,0.4)
  ) +
  theme(
    panel.grid.major = element_line("gray95", size = 0.1),
    # putting label closer to axis bc exponent makes it bigger 
    axis.title.y = element_text(margin = margin(r = -9)) 
  ) +
  xlab("time (h)") +
  ylab(expression(paste("concentration (mmol",~cm^-3,")"))) +  
  coord_cartesian(clip = "off"# to allow for plotting outside axes

panel_E

小编有话说

本文主要学到的知识点如下:

  1. 使用 geom_point() 绘制散点图, geom_ribbon()绘制丝带形状图;
  2. 使用 stat_function() 添加函数曲线;
  3. 使用 geom_errorbar()geom_errorbarh 添加误差棒(纵向与横向);
  4. 使用 coord_cartesian(clip = "off") 允许展示外轴的图形。

看完这篇文章,相信你的技能包又多了点东西。记得实操噢!看了不代表就会了~ 如果觉得内容有用的话,小编写的有心的话。给小编来杯咖啡吧!

参考资料

[1]

GitHub - marco-meer/scifig_plot_examples_R: Scientific publication figure plotting examples with R: https://github.com/marco-meer/scifig_plot_examples_R

[2]

文章: https://ggplot2.tidyverse.org/reference/theme.html


推荐: 可以保存以下照片,在b站扫该二维码,或者b站搜索【庄闪闪】观看Rmarkdown系列的视频教程。Rmarkdown视频新增两节视频(写轮眼幻灯片制作)需要视频内的文档,可在公众号回复【rmarkdown


R沟通|Rmarkdown教程(4)


R沟通|Rmarkdown教程(3)


R沟通|Rmarkdown教程(2)




您可能也对以下帖子感兴趣

文章有问题?点此查看未经处理的缓存